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Radio  frequency  based  discharges  at  atmospheric  pressures  are  the  focus  of  increased  interest  in 
aerodynamics  because  of  a  wide  range  of  potential  applications  including  specifically  actuation  in  flows  at 
moderate  speeds.  Recent  literature  describing  promising  experimental  observations,  especially  on  separation 
control,  have  spurred  efforts  in  the  development  of  parallel  theoretical  modeling  to  alleviate  limitations  in 
current  understanding  of  the  actuation  mechanism.  The  present  effort  builds  on  a  recently  developed  finite 
element-based  one-  and  two-dimensional  multi-fluid  formulation  of  plasma-sheath  for  atmospheric  optical  glow 
discharge  in  partially  ionized  gas.  The  model  was  relatively  straightforward  but  formed  the  foundation  of  a 
versatile  first-principles  based  methodology.  Higher-fidelity  models  are  included  to  yield  a  more  sophisticated 
framework  to  predict  discharge-induced  momentum  exchange.  Here,  the  complete  problem  of  a  dielectric  barrier 
discharge  based  separation  control  with  axially  displaced  electrodes  is  simulated  in  a  self-consistent  manner. 
Model  predictions  for  transient  evolution  of  charge  densities,  the  electric  field,  eleetrodynamic  force  and  induced 
gas  velocity  distributions  for  helium  gas  in  quiescent  condition  are  shown  to  mimic  trends  reported  in  the 
experimental  literature.  For  the  first  time,  resuits  also  document  the  decay  process  of  a  separation  bubble  formed 
due  to  flow  past  a  flat  plat  inclined  at  12  deg  angle  of  attack.  This  effort  sets  the  basis  for  extending  the 
formulation  further  to  include  polyphase  power  input  in  multi-dimensional  settings,  and  to  apply  the  simulation 
method  to  flows  past  common  aerodynamic  configurations. 


Nomenclature 

A,B  Coefficients 

C  Capacitance,  farad 

D  Diffusion  coefficient,  cm2/s 

d  Characteristic  length,  cm 

E  Electric  field,  V/cm 

e  Electron  charge,  coulomb 

e  Permittivity,  farad/m 

F  Solution  residual 

4>  Potential,  V 

T  Flux,  cm"2/s 

I  Current,  Amp/cm: 

/W)e  Electron  Debye  length,  '/m 

/,,  Ion  mean  free  path,  cm 

m  mass,  kg 

p  Mobility,  cm2V"1s"1 

n  Number  density,  cm"3 

to  Applied  frequency,  radians 

Q  Computational  domain 
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I.  Introduction 

ECENT  research  efforts1'7  have  demonstrated  the  feasibility  of  utilizing  dielectric  barrier  discharges  (DBDs)  to 
improve  the  lift-drag  characteristics  and  inhibit  phenomena  such  as  stall  on  airfoils.  A  particularly  attractive  feature 
of  these  plasma-based  devices  is  their  capacity  to  operate  at  atmospheric  conditions,  without  uncontrolled  macroscopic 
breakdown.  Experimental  efforts1"'  have  identified  many  of  the  key  parameters,  most  prominent  among  which  are  the 
geometric  configuration  and  the  form  of  the  applied  excitation.  Specifically,  the  embedded  electrode  is  displaced  in  the 
streamwise  direction  relative  to  the  exposed  electrode  which  lies  on  the  surface  of  the  object,  such  as  the  airfoil.  The 
asymmetry  in  the  location  of  the  electrodes,  coupled  with  the  phase  shift  of  the  electrode  when  multiple  devices  are 
present,  yields  a  directional  asymptotic  “push”  on  the  bulk  gas. 

A  typical  schematic  of  such  actuators  is  shown  in  Figure  1.  The  electrode  arrangement  is  such  that  the 
insulator  surrounds  the  grounded  electrode  and  a  voltage  powered  at  radio  frequency  (if)  is  applied  to  the  electrode 
exposed  to  the  gas.  In  another  arrangement,  both  electrodes  are  powered  separated  by  a  beat  frequency.  A  complex 
unsteady  interaction  occurs  between  the  two  electrodes,  details  of  which  depend  on  frequency,  voltage,  geometric 
configuration  and  dielectric  constants  of  the  media.  Experimental  evidence  shows  that  there  is  no  runaway  state  for  the 
parameters  under  consideration  and  that  an  asymptotic  (quasi)  periodic  state  is  reached,  with  a  dominant  frequency  that 
is  similar  to  the  input  perturbation.  For  a  given  inter-electrode  distance,  as  the  applied  voltage  becomes  sufficiently 
large,  the  dielectric  surface  adjacent  to  the  rf  electrode  produces  a  barrier  discharge,  which  weakly  ionizes  the 
surrounding  gas.  The  combination  of  electrodynamic  body  force  and  collisional  processes  whose  detailed  mechanics  is 
yet  to  be  understood,  ultimately  transfers  momentum  acquired  from  the  electric  field  by  the  charged  particles  to  the 
neutrals  which  are  the  primary  species. 
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Figure  1.  Schematic  of  the  RF  induced  dielectric  barrier  discharge  device.  The  shaded  region  is  an  insulator  with 
embedded  (grounded,  ((>  =  0)  and  exposed  (powered,  <j>0  f(cot))  electrodes  shown  in  thick  solid  line  segments. 

Benefits  of  these  actuators  include  precise  application  of  large  electrodynamic  body  forces,  active  control  of 
these  forces,  and  their  striking  effects  on  the  flow.  These  capabilities  added  with  the  absence  of  any  moving  components 
make  these  actuators  very  useful  for  flow  control  in  wall  layers  or  separated  layers  at  both  low1'7  and  high  speeds.8'9  The 
parameters  employed  in  experimental  observations  include  peak-to-peak  voltage  between  2-20  kV  at  1-50  kHz  rf,  which 
are  suitable  for  actuation  at  atmospheric  pressure  at  low  speeds  to  0(10)  torr  at  high-speeds.  Specifically  at  high 
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pressures,  the  fluid  is  highly  collisional,  inducing  an  efficient  momentum  transfer  between  charged  and  neutral  species. 
The  collisional  rf  sheath  dynamics  of  near-surface  regions  are  fundamentally  different  than  that  under  direct  current  (dc) 
and/or  collisionless  conditions  and  substantially  more  difficult  to  simulate  because  of  their  unsteady  nature.  While  de¬ 
based  methods  are  a  useful  study  for  energy  interactions  of  an  already  ionized  flow,  recent  literature8"10  shows  that 
energy  budgets  will  depend  crucially  on  dynamic  non-equilibrium  ionization  techniques.  The  load  factor  for  dc  sheath 
application  is  of  order  1,  far  from  the  most  energy  efficient  Stoletow  point,  thus  generally  unsuitable  for  ionization 
purposes.  The  experimental  successes  of  rf  induced  dielectric  barrier  discharge  (DBD)  using  air  at  atmospheric 
condition1,3  thus  ushers  tremendous  interest  in  the  electrogasdynamic  flow  control  and  energy  management  community. 
Microfilaments  of  nanosecond  duration  with  many  current  pulses  in  a  half  cycle  maintain  the  optical  glow  in  DBD.  The 
small  time  scale  of  the  discharge  limits  charge  transport  and  allows  minimal  heating  of  the  bulk  gas.6,11  Thus  most  of  the 
electric  energy  is  utilized  for  exciting  the  carrier  gas  resulting  in  efficient  applications  to  boundary  layer  flow  actuation. 

This  paper  focuses  on  modeling  the  physics  of  a  plasma  actuator,  shown  in  Figure  1,  about  which  a  surface 
barrier  discharge  controls  the  surrounding  neutral  gas  flow.  Several  configurations  have  been  proposed  and  tested  to 
exploit  the  complex  interaction  between  the  electric  field  and  the  fluid  around  such  actuators.  A  detailed  discussion  may 
be  found  in  Ref.  1,3. 

While  experimental  efforts  on  atmospheric  DBD  interactions  for  flow  actuation  are  ample  and  extensive,  the 
theoretical  modeling  effort  that  span  a  range  of  phenomenological  to  rudimentary  first  principle  based  methods2'4"712  is 
yet  to  explain  the  underlying  physics.  Roth''  postulated  a  balance  between  the  electrostatic  force  and  the  hydrodynamic 
force  to  estimate  the  induced  gas  velocity.  Enloe  et  al.2  analytically  showed  that  this  postulation  only  holds  for  a  simple 
one-dimensional  configuration  and  demonstrated  the  need  for  two  dimensional  electric  force  modeling  for  understanding 
the  DBD  flow  actuation  process.  Based  on  experimental  observation,  Shyy  et  al.4  assumed  a  triangular  distribution  of 
the  electrostatic  body  force  downstream  of  the  electrode  with  a  constant  charge  density  and  then  utilized  the  Navier- 
Stokes  equation  to  predict  its  effects  on  the  flow.  Recently,  Gaitonde  et  al.12  have  described  the  response  of  the  flow  past 
a  stalled  NACA  0015  airfoil  at  15  deg  angle  of  attack  to  phenomenological  approximate  body  forces  originating  from 
asymmetric  dielectric-barrier-discharge  actuators  using  direct  numerical  simulations. 

For  accuracy  and  fidelity,  however,  it  is  imperative  that  the  force  model  be  derived  from  first  principles  through 
a  simulation  of  the  elementary  mechanisms  that  yield  the  discharge  characteristics  and  flow  actuation  similar  to  that 
found  in  the  experiments.  Roy  and  Gaitonde  demonstrated  such  a  model,  first  for  volume  discharge  between  two 
dielectric  coated  electrodes5  and  then  for  capturing  effects  of  surface  discharge  on  quiescent  helium  gas6’7,  using  a  finite 
element-based  multi-dimensional  multi-fluid  formulation  of  plasma-sheath  at  atmospheric  conditions.  While  the  first 
model  was  relatively  basic,  the  two  dimensional  model  of  collisional  surface  dielectric  barrier  discharge  documented  a 
consistent  first  principles  formulation  of  charge  and  neutral  number  densities,  their  momentum  dynamics,  electrostatic 
field  and  potential  distribution  for  an  asymmetric  DBD  configuration.  The  model  treated  the  insulator  and  the  gas 
simultaneously  and  integrated  the  Poisson  equation  into  the  charge  and  gas  dynamics.  The  study  was  focused  on 
understanding  the  effects  of  symmetric  and  asymmetric  electrode  arrangements.  Simulation  results  for  asymmetric 
arrangement  from  Ref.  6  as  shown  in  Figure  2  compared  qualitatively  with  reported  rf  discharges1"3  in  partially  ionized 
helium  gas. 


Figure  2.  Plasma  actuation  of  quiescent  helium  gas.6  (a)  Charge  particle  distribution,  (b)  force  variation  about  the 


electrode-dielectric  surface,  (c)  streamwise  gas  velocity  profiles  at  different  locations  along  the  flow. 

As  a  logical  extension  of  our  previous  work,  the  goal  of  the  present  effort  is  to  implement  the  developed  two- 
dimensional  model  to  practical  single  electrode-pair  systems  for  understanding  the  response  of  body  force  characteristics 
resulting  from  the  collisional  surface  discharge  to  separated  gas  flows  at  atmospheric  pressures.  An  efficient,  robust, 
module-based  multiscale  ionized  gas  (MIG)  flow  finite-element  code  suitable  for  accurate  imposition  of  complicated 
boundary  conditions  is  employed.  The  body  force  is  computed  self-consistently  from  the  three-body  electro  gas 
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dynamics  using  the  applied  rf  potential  and  the  dielectric  properties  of  insulator  and  gas  media  based  on  first  principles. 
For  simplicity,  here  again  helium  is  used  as  the  carrier  gas  because  of  its  relative  stability  in  maintaining  the  discharge 
optical  glow  (as  compared  to  air,  nitrogen,  oxygen  and  argon  which  transition  into  an  unstable  filamentary  discharge  for 
the  same  electrode  gap  and  pressure").  Note  the  model  is  easily  extendible  to  complex  configurations  and  more  accurate 
transport  properties  when  they  become  available  from  ongoing  experimental  efforts  being  performed  elsewhere.  Two 
simulations,  one  with  quiescent  gas,  the  other  with  an  imposed  upstream  gas  velocity  are  documented  here  showing  not 
only  plasma  events  but  also  the  evolution  of  separated  flow  control  about  an  actuator. 

This  paper  is  organized  as  follows.  Section  II  describes  the  problem  statement.  Section  III  specifies  the 
boundary  conditions.  Section  IV  explains  the  methodology  used  to  solve  the  system  of  equations.  Section  V  describes 
and  interprets  the  computed  results.  Section  VI  summarizes  the  conclusions. 


II.  Problem  Statement 


Assuming  only  positive  ions,  negative  electrons  and  neutral  particles,  the  following  two-dimensional  three- 
body  collisional  plasma-sheath  model  is  solved. 


Charge  continuity: 

dna  dnaV^ 

— —  H - -  =  n  z  —  rn  n . 

(la) 

dt 

e  e  i 

dXj 

Charge  momentum: 

n 

v .  = 

a  eg 

-sgn  (e)nana^~Dad^- 

OXj  OX  j 

(lb) 

Potential: 

£ 

dx) 

=  e(ne-ni) 

(2) 

V  J  J 

Neutral  continuity: 

dn„  ,  dn«vnj 

— -  + - -  =  —n.z  +  rn 

(3) 

dt 

dXj 

The  charge  particle  a  =  e,i  distributions  are  considered  non-Maxwellian.  The  electron  temperature  is  nearly 
uniform  at  leV  =  1 1,600K  and  the  ions  and  neutrals  are  in  local  thermal  equilibrium  at  300K.  The  working  gas  helium  is 
maintained  at  300  torr  bulk  pressure  at  300K.  For  the  operating  condition  of  interest,  the  inertia  terms  are  small  relative 
to  the  effect  of  the  electrostatic  field  and  collisional  interactions,  and  thus  neglected.  The  collisional  momentum 
exchange  is  vital  to  the  charge-to-neutral  momentum  transfer. 

The  electron  diffusion  is  obtained  from  Einstein  relation,  £)  =  {Te  I e) jue,  where  Te  is  the  energy  in  electron 

volts,  e  is  the  elementary  charge,  s  is  the  permittivity,  and  / ue=e/(meve/, )  is  mobility  of  an  electron,  where  ve„*40l::/s  is 
the  electron-neutral  collision  frequency.  In  this  context,  we  note  that  transport  properties,  which  have  been  taken  from 
the  literature,  are  to  be  viewed  as  nominal  values  facilitating  the  development  of  the  numerical  framework,  rather  than  to 
fill  gaps  in  current  thermo-chemical  data  for  the  environment  of  interest.  In  this  same  vein,  D,  is  the  ion  diffusion14  in 
cm2/s  at  300K,  and  the  ion  mobility  //,  is  given  as  an  empirical  relation  that  contains  the  electric  field  and  pressure, 
assuring  consideration  of  collisional  effects.7 14  Flowever,  quantitative  results  will  depend  on  the  nominally  chosen 
values.  The  ionization  rate  z  for  the  working  gas  is  z‘/nn  =  <Ved(V^>  ,  the  averaging  is  done  over  the  velocities  of  the 
electrons  whose  energy  is  sufficient  for  ionization  mVe2/2  >//„  the  first  ionization  potential.  Flere  z  is  based  on  the 
Townsend  first  ionization.6'7  The  ionization  process  can  be  described  as  e  +  He  — >  He+  +  e  +  e’  where  e  and  e’  have 
different  energy  levels.  Processes  like  He—>He++  and  He'  —>He++  may  also  play  an  important  role  for  high  voltages,  but 
are  not  considered  here.  The  coefficient  of  recombination  is  given  as,  r  =  lv elh<7ra(V elh  )\  =  1.09  x  lO~uT~9l2ne  cm3/s.  Flere, 


V,ih  is  the  electron  thermal  velocity  and  aei  is  the  electron-ion  collision  cross-section.  For  simplicity,  we  have  assumed 
no  secondary  electron  emission. 

Flere,  the  electron  plasma  frequency,  ft^,e  =  2e^l(Kne/me),  is  much  higher  than  the  applied  voltage  frequency. 
Also,  the  momentum  exchange  terms  due  to  electron-electron  and  ion-ion  interaction  are  negligible  in  comparison  with 
the  electron-ion  and  ion-neutral  momentum  exchange  as  the  relative  drift  between  similar  particles  is  small  in 
comparison  with  the  drift  between  electrons,  ions  and  neutrals. 

The  dynamics  of  neutral  particles  are  determined  using  the  following  neutral  momentum  equations: 


dV  ■ 

w 


+  v. 


8V  ■ 

»J  _ 


m  n , 


( 1 v;  - v„ )  +  zrr' ( v,  -  r, ) + v v,  /"• 


(4) 


dt  "  dx,  m  n ,  '  "  ' '  m  n 

Inn  n 

The  effect  of  pressure  on  quiescent  neutral  flow  is  assumed  negligible.  This  was  chosen  for  separating  the 
effect  of  electric  field  from  the  pressure  field.  This  assumption  is  also  reasonable  for  high  Reynolds  number  flows.  At 
this  stage,  any  balance  between  the  electrostatic  and  fluid  dynamic  pressure  gradients  cannot  be  guaranteed  for  a 
fundamentally  unsteady  process.  We  use  the  electrodynamic  body  force  e(/7,-”e)E  (although  no  Lorentz  force  is  present. 
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the  force  is  created  by  charge  transport)  to  compute  the  averaged  velocity  of  the  fluid  is  based  on  the  following 
continuity  and  momentum  equations: 


dp /dt  +  dpV f]  jdxj  =  0 


(5a) 


dVr  dV„  uj  ,  s 
—  +  V„^  =  -e(ni-nc,)Ei 
dt  fl  dx,  p  1 


R  T  dp  d 

- 1 - T  7  > 

p  dxj  dx,  ' 


(5b) 


where  T  =tL  +  with  pas  the  gas  viscosity  (=2.066xl0-5  Ns/m2)  and  5,7  is  the 

Jl  P  \  Sxj  dx,  J  3  p  8xk  Jl 


Kronecker  delta.  The  factor  tn  is  introduced  to  modulate  the  effect  of  electric  field.  Helium  is  assumed  as  an  ideal  gas. 

The  multiscale  ionized  gas  (MIG)  flow  code  was  used  to  solve  Eqns.  (l)-(5).  MIG  has  a  module  based 
computational  platform  that  implements  the  high-fidelity  finite-element  (FE)  procedure  adapted  from  fluid  dynamics  to 
overcome  the  stiffness  of  the  equations  generated  by  multi-species  charge  separation  phenomena.  This  code  has  been 
used  and  validated  for  several  fluid  dynamics15,  ionized  gas5"7’1617  and  micro/nanoscale ls' 19  flow  problems.  FE 
techniques  are  especially  suitable  for  their  adaptability  to  arbitrary  multidimensional  geometries  and  boundary 
conditions.  Here,  a  2D  biquadratic  FE  formulation  is  employed  along  with  4th  order  Runge-Kutta  time  marching  to  solve 
the  continuity  and  momentum  equations  for  all  species  and  fluid.  The  solution  process  consists  of  two  steps.  First,  using 
Eqns.  (la),  (2),  (3)  and  (4),  a  global  matrix  is  formed  and  solved  simultaneously  obviating  the  need  for  any  special  sub¬ 
iteration  for  the  Poisson  solver.  The  species  density  and  charge  velocity  thus  calculated  are  then  used  in  Eqn.  (5a-b)  to 
predict  the  fluid  density  and  velocity.  The  Galerkin  weak  statement  associated  with  a  variational  integral  underlines  the 
development  of  this  numerical  algorithm.  Note  that  the  selection  of  test  function  orthogonal  to  the  trial  function  in 
Galerkin  formulations  guarantees  the  minimization  of  error20,  making  it  better  suited  for  these  problems.  Details  of  the 
code  are  reported  elsewhere  (for  example,  see  Ref.  15,16).  For  the  species  equations  (l)-(4),  a  solution  adaptive 
timescale  ofupto  100  times  the  dielectric  relaxation  timescale21  is  adopted  to  calculate  the  charge  separation  and  electric 
field  up  to  a  periodic  asymptote.  The  gas  equations  (5a-b)  are  then  solved  using  the  volume  specific  electrodynamic 
body  force  to  predict  the  fluid  velocity  and  continuity  evolving  in  time  with  the  flow  time  scale  of  ~ms.  The  solution  of 
the  Newton-Raphson  iteration  is  converged  at  any  given  timestep  when  the  maximum  value  of  the  residual,  relative  L2  - 
norm  for  each  of  the  state  variables,  becomes  smaller  than  a  chosen  convergence  criterion  of  10'4. 


III.  Boundary  Conditions 

The  two  dimensional  computational  domain  (-1.252,  3.880cm)x(-0.1,  2.0cm)  consists  of  a  Kapton  insulator  in 
the  lower  half  (y:  -0.1,  0cm)  with  a  dielectric  constant  £^=3.5  eq  and  the  upper  half  (y:  0,  2cm)  filled  with  inert  helium 
gas  of  £f=  1.0055  s0,  where  s0  is  permittivity  of  vacuum.  Inside  the  insulator  the  current  due  to  motion  of  charged 
particles  is  forced  to  zero  while  the  displacement  current  is  balanced  with  the  total  current  at  the  gas-dielectric  boundary. 
The  schematic  in  Figure  1  shows  two  electrodes  in  which  the  bottom  electrode  is  grounded  and  an  rf  alternating 
frequency  of  5  kHz  with  rms  potential  of  1  kV  is  imposed  at  the  top  electrode.  Each  electrode  is  infinitesimally  thin  and 
1.2  cm  long.  They  are  vertically  displaced  by  1mm  with  no  horizontal  overlap.  The  highly  stretched  computational  mesh 
shown  in  Figure  3  consists  of  51  x  53  two-dimensional  bi-quadratic  Lagrange  finite  elements.  The  insulator  material  is 
modeled  using  12  x  53  biquadratic  elements.  The  nodes  are  attracted  towards  the  wall  and  electrode  edges  for  capturing 
better  solution  details.  The  highest  and  lowest  elemental  aspect  ratios  are  approximately  90  and  10,  respectively. 


Figure  3.  Computational  mesh  shows  highly  non  uniform  nature. 
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In  order  to  identify  the  charge  characteristics  and  its  effect  on  gas  flow,  two  different  simulations  were 
performed  on  a  single  geometric  configuration  with  an  asymmetrically  displaced  electrode-pair  arrangement.  In  the  first, 
no  gas  velocity  was  imposed,  i.e.,  quiescent  flow.  Hereafter,  we  will  denote  this  as  Case  A.  In  the  second,  an  upstream 
flow  of  12  m/s  was  imposed  at  x  =  -1.252  cm  with  an  angle  of  attack  (AOA)  of  +12°.  This  will  be  noted  as  Case  B.  For 
all  simulations  electrons  are  assumed  to  be  isothermal  at  the  boundary  and  maintained  at  leV  (~1 1,600K)  while  the  ions 
are  cold  (300K)  at  300  torr.  Zero  flux,  i.e.,  homogeneous  Neumann  boundary  conditions  are  applied  at  all  other 
boundaries.  The  solutions  are  verified  by  qualitative  comparison  with  the  reported  results.  The  results  are  then  employed 
to  explore  the  enhancement  of  near  wall  fluid  velocity.  No-slip  gas  velocity  boundary  condition  is  imposed  for  all  cases. 
Also  the  insulator  is  considered  impermeable  to  gas. 


IV.  Methodology 

The  ionized  helium  gas  is  numerically  modeled  using  the  finite-element  based  MIG  flow  code.  The  code  is  modular  and 
separate  subroutines  can  be  written  to  model  different  physics.  Here,  the  equation  sets  (1)  -  (5)  can  be  written  with 
operator  L  as  L  (q)  =  0  ;  where  q  contains  state  variables  like  densities,  velocities  and  potential.  Multiplying  with  a 
permissible  test  function  77  and  integrating  over  the  spatially  discretized  domain  Q,  the  variational  statement  results  in 
the  weak  form20 


WS  =3. 


|  [/jL(q)dr] 


=  0 


for  a  discretization  h  of  Q  =  1JQ(  and  3e  is  the  non-overlapping  sum  over  the  elements.  Thus,  for  example,  the  GWS 
form  of  Eq.  (2)  becomes, 

X  j  +  j  Wdz{nt}t  -  j  >l>lTdz{nj}e  -  J  77 ^-dz{<t>}e  =  F$ 


dz 


(6) 


where  is  the  solution  residual. 

The  Jacobian  matrix  J=\dF/dq\  in  the  global  \J\.{8q}  =  -\F\  is  resolved  using  Generalized  Minimal  RESidual 


solver  with  LU-preconditioner  for  updating  change  in  discretized  solution  vector  q  at  each  iteration.  The  Newton- 
Raphson  scheme  is  used  to  formulate  the  matrix  algebra  resulting  from  nonlinear  ordinary  differential  equation  (ODE) 
systems.  The  GMRES  solver  is  well  suited  to  handle  the  sparseness  of  the  resulting  stiff  matrix.  The  solution  is  assumed 
to  have  converged  at  a  time-step  when  the  L2  norm  of  all  solution  variables  and  residual  are  below  a  chosen  convergence 
criterion  (£•).  The  convergence  criterion  for  all  variables  at  any  iteration  is  set  at  10"4.  Solution  stability  is  ensured  by 
appropriate  selection  of  time  marching  step  size  and  the  introduction  of  artificial  diffusion. 


V.  Results  and  Discussion 

Figures  4  through  11  present  computed  results  for  charge  densities,  electric  field  lines,  current  and  gas 
dynamics  for  simulated  cases.  Figure  4  displays  unsteady  evolution  of  predicted  force  vectors  with  electric  field  lines 
overlaid  on  charge  density  contours  (cm'3)  for  the  quiescent  Case  A.  The  electric  field  lines  trace  the  vector  path  moving 
in  a  trajectory  from  the  exposed  instantaneous  anode  to  the  grounded  cathode  inside  the  dielectric  in  an  asymmetric 
fashion  showing  a  directional  bias,  first  (at  jt/2  radians,  Fig.  4a)  towards  the  right,  then  towards  the  left  (at  n,  Fig.  4b  and 
at  3n/2  radians,  Fig.  4c),  and  finally  towards  the  right  again  at  2 71  radians  (Fig.  4d).  The  computed  results  are  similar  to 
the  experimental  data  showing  that  the  exposed  electrode  is  situated  upstream  of  the  peak  location  of  the  electric  field 
and  thus  the  imparted  momentum.  The  dominant  acceleration  of  ions  following  these  lines  induces  the  neutral  gas 
particles  through  a  mechanism  entwined  with  the  body  force  and  the  ion-neutral  collisional  momentum  transfer.  In  Fig. 
4a,  the  charge  density  contours  shows  a  positive  peak  of  3.4xl09  cm"3  above  the  embedded  electrode  (a  few  Debye 
length  away  from  the  surface)  and  a  negative  peak  of  -1.9xl09  cm'3  on  the  exposed  electrode  showing  the  effects  of 
electric  field.  The  ion  density  dominates  the  actuator  surface  at  n  radians  with  a  peak  of  6.3xl09  cm"3.  At  3n/2,  the 
negative  charge  (-2.7xl09cm"3)  accumulates  downstream  of  the  exposed  electrode  and  charge  density  extrema  gradually 
decay  as  the  cycle  moves  to  2n  radians. 

The  body  force  calculated  based  on  the  product  of  charge  and  electric  field  shows  very  strong  forward  and 
downward  moving  vectors  downstream  of  the  exposed  electrode  for  the  positive  part  of  the  cycle  (see  Fig.  4a).  At  the 
negative  peak  of  the  cycle  in  Fig.  4c,  the  force  vectors  concentrate  just  upstream  of  the  right  edge  of  the  powered 
electrode  with  majority  forces  going  backward  and  downward.  However,  for  the  rest  of  the  cycle  the  force  is  much 
smaller  as  compared  to  the  positive  part  (-10%,  42%  and  2%  of  the  peak  magnitude  respectively  at 
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ji,  3ji/2  and  2n  radians).  The  average  body  force  will  thus  push  the  bulk  gas  in  the  forward  and  downward  direction. 
Recent  numerical  simulation6'7  documented  the  gas  velocity  downstream  of  the  right  edge  of  the  exposed  anode  with  a 
two  orders  of  magnitude  smaller  negative  velocity  near  the  left  edge  of  the  exposed  electrode  similar  to  what  we  observe 
here. 


(d)  2 n  rad 

Figure  4.  Time  evolution  of  computed  electric  field  lines  overlaid  on  charge  density  contours  (cm-3)  with  near 

electrode  volume  specific  force  details. 
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Figure  5  plots  the  streamwise  component  of  the  time  average  of  volume  specific  body  force  for  Case  A.  The 
line  trace  of  the  force  vectors  is  showing  a  directional  bifurcation  (similar  to  that  reported  in  a  recent  PIV  experiment22, 
Fig.  2)  just  downstream  of  the  exposed  electrode.  This  will  tend  to  guide  the  flow  in  both  direction  just  at  about  x  =  2cm 
location.  The  positive  peak  very  close  to  the  dielectric  surface  is  three  times  higher  than  the  negative  peak.  Due  to  fluid 
inertia,  the  bulk  gas  will  only  respond  to  this  average  force  which  will  ensure  its  net  forward  motion.  The  momentum 
thus  imparted  to  the  gas  will  induce  a  velocity  along  the  dielectric  surface.  Figure  6  plots  the  streamwise  component  of 
the  gas  velocity  computed  from  Eq.  (5)  for  Case  A  at  six  local  vertical  line  plots  upstream  and  downstream  of  the 
electrode  edge  and  shows  a  wall  jet  like  feature.  The  zero  flow  initial  condition  makes  the  computational  problem  more 
challenging.  The  peak  of  this  wall  jet  of  2.4  m/s  occurs  at  6  mm  downstream  of  the  exposed  electrode  beyond  which 
diffusive  momentum  transfer  increases  the  height  of  the  region  influenced  by  the  DBD,  and  the  peak  velocity 
diminishes.  Similar  velocity  profdes  have  been  observed  by  several  experimental  groups.1'3'22  Note  that  for  helium  gas 
the  mobility  is  much  higher  than  air  and  hence  the  induced  velocity  is  higher. 

The  peak  electron  current  (not  shown)  is  observed  at  the  two  corners  (right  and  left)  of  the  exposed  electrode 
with  a  peak  value  of  (910 2  mA  at  n  radians  and  the  minimum  current  (two  orders  of  magnitude  less  than  the  peak)  is 
computed  at  jr/2  radians  of  the  phase.  An  explanation  for  this  behavior  may  be  proposed  based  on  arguments  pertaining 
to  the  observed  self-limiting  nature  of  the  discharge.  In  surface  DBD  there  is  no  real  current  path  between  the  electrodes. 
Flowever,  the  charge  build  up  on  the  insulator  surface  opposes  the  applied  voltage,  shutting  of  the  discharge  just  beyond 
the  peak  when  the  applied  voltage  starts  decreasing,  thus  self-limiting  the  system.  This  is  also  confirmed  by  the 
experimental  data  (see  for  example  Fig.  2  of  Ref.  23)  where  the  current  magnitude  peaks  at  about  n  radians. 


Figure  5.  Time  average  of  streamwise  component  of  the  pN  force  per  unit  volume  about  the  surface  of  the 
actuator  shows  the  dominance  of  the  streamwise  forward  (positive)  force  component.  Line  traces  show  the 

direction  of  force  which  will  guide  the  local  flow. 


Figure  6.  Computed  streamwise  velocity  induced  in  a  quiescent  helium  gas. 

For  the  actuator  model  to  work  for  design  purposes,  it  is  essential  to  test  it  for  low  speed  flow  conditions 
where  reasonable  amount  of  supporting  experiments  have  been  reported.  Figure  7  shows  the  transient  effect  of  DBD 
actuator  for  realistic  gas  flow  control.  The  timescale  of  gas  flow  for  the  present  dimensions  is  in  milliseconds  while  the 
plasma  timescale  is  several  orders  of  magnitude  less.  For  this  Case  B  simulation,  first  an  initial  flow  condition  is 
generated  over  the  flat  plate  kept  at  a  12°  angle  with  respect  to  the  helium  gas  inflow  without  the  rf  power  turned  on. 
Figure  7a  plots  the  stream  traces  of  gas  particles  based  on  the  velocity  vectors  of  the  initial  flow  field  and  shows  a  strong 
separated  flow  (see  inset)  just  downstream  of  the  leading  edge  of  the  plate.  Then  the  power  is  turned  on  for  the  same 
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electrode  arrangement  as  presented  earlier.  Figures  7b, 7c  and  7d  shows  the  gradual  removal  of  the  recirculation  bubble 
due  to  the  separation  as  a  function  of  time.  The  strong  positive  force  contours  (in  red)  shows  the  time  average  dynamics 
of  streamwise  component  as  it  literally  pushes  the  recirculation  off  the  plate  within  6ms.  The  fluid  reacts  to  the 
generated  mean  body  force  but  there  may  be  a  slow  interaction  with  the  local  plasma  generation  as  time  progresses.  To 
our  knowledge,  such  description  has  not  been  shown  in  the  reported  literature  to-date.  Specific  details  of  gas  velocity 
components  and  the  streamwise  force  is  presented  in  Figure  8  at  five  locations  along  the  dielectric  at  Z=6ms.  It  is 
interesting  to  note  that  the  peak  force  is  located  at  x  =  0.7  cm  in  the  middle  of  the  initial  bubble  while  the  peak  of  the 
wall  jet  is  occurring  downstream  of  the  exposed  electrode  at  x  =  2.0  cm.  The  location  of  the  actuator  is  thus  important  to 
control  the  flow  effectively. 


O 


O 


0  2 
Streamwise  (cm) 

(a)  t  =  Os 


o 
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O 
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0  2 
Streamwise  (cm) 

(b)  t  =  2  ms 


0  2 

Streamwise  (cm)  Streamwise  (cm) 

(c)  t  =  5  ms  (d)  t  =  6  ms 

Figure  7.  Separation  control  for  an  incoming  flow  with  +12°  AOA  over  a  flat  plate  shows  the  gradual  flow 
attachment  eradicating  the  separation  bubble  in  a  few  milliseconds.  The  velocity  streamtraces  are  plotted  over 
the  streamwise  specific  force  contours  showing  a  strong  average  force  over  the  exposed  electrode. 


Figure  8.  Velocity  components  and  streamwise  volume  specific  force  distribution  at  specific  locations  (x  in  cm) 

normal  to  the  dielectric  surface  after  6  ms. 
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VI.  Conclusions 

A  two-dimensional  finite  element  based  formulation  of  plasma-fluid  interactions  is  given  for  partially  ionized  plasma 
using  the  multi-component  fluid  equations.  The  model  is  applied  to  simulate  an  atmospheric  surface  dielectric  barrier 
discharge  for  partially  ionized  helium  gas.  The  computed  results  are  similar  to  the  experimental  data  showing  that  the 
exposed  electrode  is  situated  upstream  of  the  peak  location  of  the  electric  field  and  thus  the  imparted  momentum.  The 
two-dimensional  predictions  mimic  the  self-limiting  discharge  mechanism  due  to  the  charge  accumulation  on  the 
insulator  surface  and  the  lowest  electron  current  just  beyond  the  peak  applied  voltage.  Simulation  results  also  show  the 
effectiveness  of  plasma  actuators  for  control  of  largely  separated  flows  useful  for  many  practical  applications  including 
airfoils  and  turbine  blades.  More  investigation  is  needed  to  ascertain  optimum  electrochemical  parameters.  In  our 
simulations,  the  electric  field  varies,  but  reaches  a  quasi-periodic  asymptote.  The  electrode  placement  plays  a  significant 
role  in  controlling  the  separated  flow.  We  believe  the  fluid  reacts  to  the  generated  mean  body  force  but  may  also  weakly 
interacts  as  time  evolves,  however,  this  comment  requires  further  clarification.  In  the  near  future,  the  model  will  be 
extended  to  air  with  negative  ions  and  additional  new  mechanisms  in  the  source  terms.  A  model  for  realistic  fully  three- 
dimensional  geometric  and  electrode  configuration  is  also  under  development,  as  is  an  exploration  of  the  effect  of 
different  voltage  shapes  (e.g.  saw  tooth,  square  wave  etc.).  The  present  effort  thus  provides  a  practical  tool  to  augment 
experimental  observations  in  exploring  flow  control  concepts  and  in  developing  suitable  inputs  for  traditional  fluid 
dynamics  codes  based  on  the  Navier-Stokes  equations. 
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